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I. INTRODUCTION 

While providing a formally simple solution for the 
quantum dynamics of a physical system, the Feynman 
path integral methodi generates one of the most dim- 
cult problems in computational physics, when it comes 
to the actual simulation on a "classical" computer: the 
dynamical sign problem,^ The highly oscillatory inte- 
grals appearing in the Feynman path integral expression 
of the propagator cannot be computed by direct Monte 
Carlo techniques, for there is no suitable importance 
function that would transform the propagator into an 
integral against a probability distribution^ Despite this 
inherent limitation, the Monte Carlo methods have been 
applied to the finite-temperature dynamics, through two 
different approaches, mainly. Both techniques alleviate 
only partly the exponential loss of signal that is indica- 
tive of the dynamical sign problem. The first approach 
tries to construct an appropriate importance function by 
convoluting the highly oscillatory integrand with a local 
distribution probability, whether a continuous^LS or a 
discrete oneiS While the research continues, actual appli- 
cations of such techniques to realistic physical systems 
are, at present, rare. 

The second approach, of which the present develop- 
ment is part, attempts to reconstruct certain dynami- 
cal correlation functions from the imaginary-time coun- 
terparts. While analytical continuation arguments show 
that this is possible in principle ; 10 : 11 the resulting algo- 
rithms always involve the resolution of certain ill-posed 
numerical problems, as for instance inverse real Laplace 
transforms or inverse moment problems. Such inverse 
problems are highly unstable and suffer from an exponen- 
tial amplification of the errors in the input data. They 
require very accurate Monte Carlo data, a careful choice 
of the type of data that is computed, as well as appro- 



priate choices of inversion and regularization algorithms. 

The most common strategy for performing the ana- 
lytic continuation is based on the inversion of a two-sided 
real Laplace transform with noisy input data, to recover 
a spectral functioniiSiiiiii^ As already mentioned be- 
fore, the inversion problem is ill-posed and, as a conse- 
quence, the inversion algorithms are highly unstable. The 
lack of continuity of the inverse Laplace transform with 
the input data causes any inversion algorithm to amplify 
the errors in the input data in an exponential fashion. 
Because these errors are of statistical nature, one is in- 
clined to believe that it is virtually impossible to recover 
any useful dynamical information. However, most no- 
tably by use of methods of Bayesian statistical inference 
with entropic priors various research groups have 
been successful in obtaining limited but useful and some- 
times surprisingly reliable quantum dynamical informa- 
tion, whether in the form of spectrarii£*i2*22i2I quantum 
rates of reaction ^22*2^4 or diffusion constants^ 

The Monte Carlo data computed for the methods 
based on the inverse Laplace transform are usually the 
values on a grid of some imaginary-time correlation 
function^ values that have proportional statistical er- 
rors and are highly redundant. However, to a larger ex- 
tent, the quality of the reconstructed spectral density is 
controlled by the errors of the relative differences, in ad- 
dition to the errors of the absolute values of these data. 
It is then apparent that ensuring low relative errors for 
such differences is a definite way of improving the quality 
of the final results, as well as the stability of the inver- 
sion algorithms. Depending on the order k of the finite- 
difference scheme considered, the value of such a differ- 
ence decreases as a polynomial of order k with the mesh 
of the grid, and so must decrease its error. It is then ap- 
parent that good quality of the input data requires good 
error bars, not only for the imaginary-time correlation 
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functions, but also for their high-order derivatives. At 
extreme, one may consider that the input data consist 
of the value of the imaginary-time correlation function 
at time zero only, together with the sequence of deriva- 
tives at origin. If such data are computed, the recon- 
struction of the spectral density involves the resolution 
of an inverse moment problem, as we shall discuss in 
the next section. The values of the derivatives at origin 
of the imaginary-time correlation function become the 
moments of the spectral density, except perhaps for a 
normalization coefficient. The inverse moment problem 
is expected to be more stable than the inverse Laplace 
transform, with respect to the relative errors in the input 
data. 

It is quite unfortunate that, for general quantum prob- 
lems, especially those in continuum space, the computa- 
tion of derivatives at origin of imaginary-time correlation 
functions is extremely difficult. When available, the mo- 
ment information can act as a stabilizing factor for the 
techniques based on the inverse Laplace transform. For 
instance, Whit o 2 ^' 2 ? utilized the first two moments of 
the spectral function for the two-dimensional Hubbard 
model as additional constraints in a maximum entropy 
approach, with remarkable success. For Fermionic sys- 
tems, Caffarel and Ceperlej* 2 ^ utilized the first moment 
(average energy) of the spectral overlap function, mo- 
ment evaluated by quantum Monte Carlo, to stabilize 
their maximum entropy computations. Another way to 
make use of a limited number of low-order moments is to 
incorporate the information into the default model. This 
technique was utilized by Diesz and coworkers to recon- 
struct spectral weight functions for the one-dimensional 
t-jS2 and Heisenberg2£ models. 

The limited use of moment information in the exam- 
ples mentioned in the preceding paragraph is due to the 
difficulties encountered in the actual computation of the 
moments: for example, only low-order moments can be 
computed in an efficient way by means of sum rules. In 
principle, when moments can be computed effectively, 
full-fledged moment techniques may be developed. How- 
ever, as far as the present author is aware, this is only the 
case for the computation of moments for certain sparse 
Hamiltonian matrices, moments that can be computed 
in O(N) operations by stochastic methods, as shown by 
Skilling^L as well as by Silver and Roderi^ 2 " The avail- 
ability of such information has led various groups to the 
application of Bayesian inference methods? 3 ^ kernel poly- 
nomial methods f 3 ^ 4 * 3 ^ or both2& to the development of 
linear scaling algorithms for the resolution of densities of 
states, in electronic structure calculations. 

There are a couple of problems in continuum space that 
would particularly benefit from a moment approach. For 
these problems, the spectral function can be made rather 
featureless and the number of necessary moments can be 
made rather small by the variation of certain physical 
parameters. One interesting case is the aforementioned 
problem of computing the Fermion ground state by quan- 
tum Monte Carlo?^ where the complexity of the spec- 



tral overlap can be greatly reduced by a proper choice 
of the antisymmetric trial function (this spectral density 
becomes a single delta function, in the limit that the ex- 
act antisymmetric trial function is utilized). In chemical 
physics, a very important problem is the computation of 
the quantum rate of reaction by time-integrating the flux- 
flux correlation function associated with a surface that di- 
vides the reactants from the productsiSSi2LS It is known 
that the quantum rate of reaction does not depend upon 
the specific choice of dividing surface, although the com- 
plexity of the flux- flux correlation function is strongly de- 
pendent upon such a choice^ Thus, provided that "the 
right" choice of surface is made, the spectral density of 
the flux autocorrelation function can be recovered effec- 
tively from a few moments, at least in principle. The 
only requirement is that these moments, or, equivalently, 
the derivatives at origin of the imaginary-time flux-flux 
correlation function, be computed with sufficient preci- 
sion. A path-integral technique that shall be presented 
in Section III is advocated as an effective way to compute 
derivatives at origin of correlation functions, for the type 
of problems discussed in the present paragraph. 

The first part of the paper provides a formal proof 
that the sequence of moments uniquely determines the 
spectral density and, therefore, the autocorrelation func- 
tion. In addition, a convergence result is proved in 
order to identify a class of reconstruction algorithms 
for the autocorrelation functions. This result, which 
is the statement of Th. [21 demonstrates that all algo- 
rithms that preserve the positivity of the underlying spec- 
tral function (such algorithms are said to be positivity- 
preserving)^ and which exactly match the first n mo- 
ments, lead to the correct autocorrelation function, in 
the limit n — > oo. The algorithms based on the maxi- 
mum entropy prmcipleiSi^SLal as well as those based on 
kernel density functions 3 ^* 3 ^ 3 ^ are examples of such algo- 
rithms. Although the proofs are conducted for thermally- 
symmetrized autocorrelation functions, Th. [3] applies for 
all correlation functions, the power spectra of which are 
positive distributions. 

The larger part of the paper is concerned with the com- 
putation of derivatives at origin of correlation functions 
for problems in continuum space. A general strategy for 
developing estimators having finite variance in the limit 
of an infinite number of path variables is discussed and 
illustrated for the case of the flux autocorrelation func- 
tion. This strategy follows the general guidance of Pre- 
descu and Doll of burying the time dependence of paths 
into the potential part of the Feynman-Kac formula. 42 
As implemented in the present paper, the computation of 
the derivatives requires the utilization of finite-difference 
schemes. Such an approach has been successfully utilized 
in recent work for the numerical evaluation of several 
thermodynamic energy and heat-capacity estimators. 43 
It is imperative to mention that no differentiation of 
Monte Carlo data is ever attempted. Rather, the finite- 
difference scheme replaces the analytical evaluation of the 
derivatives of a deterministic function, evaluation that 
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leads to expressions involving a large number of high- 
order partial derivatives of the potential, if performed. 
The additional problem we must face in the present pa- 
per is that the utilization of finite-difference schemes for 
derivatives beyond a certain order requires extended pre- 
cision arithmetic, which may be a serious programming 
nuisance. Alternative techniques, as for instance Lyness' 
method^ are possible, but they require analytic contin- 
uation of the potential in e?-dimensional complex spaces. 
Such alternatives will be investigated in future work. 

The present paper is limited to demonstrating that the 
advocated technique actually works. Any issues of effi- 
ciency are postponed for future studies. In particular, 
these studies will have to address the scaling of the vari- 
ance of the estimators utilized for the computation of 
the moments with the order of the derivatives, the di- 
mensionality of the system, and the inverse temperature. 
However, the numerical results presented in Section IV 
(these results are quantum rates for a symmetric Eckart 
barrier) show that the technique discussed in the present 
paper is a useful tool for obtaining quantum information 
of known computational difficulty. 



II. THE RECONSTRUCTION OF 
AUTOCORRELATION FUNCTIONS AS AN 
INVERSE MOMENT PROBLEM 

In physics, the quantum dynamical information mea- 
sured in experiments can generally be expressed in terms 
of quantum correlation functions of the type 
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whenever the linear response theory provides a good ap- 
proximation of the measuring physical process^ The op- 
erator H stands for the Hamiltonian of the system, a self- 
adjoint and bounded from below operator, whereas f€t 
and (3 — l/(fcgT) > are the real time and the inverse 
temperature, respectively. denotes the adjoint of the 
operator O. 

The normalization term tr (e~P H ) in Eq. is not 
relevant for our discussion and we drop it from now on. 
Using trace invariance in Eq. IjTJl. we obtain 



Co(t)=tr 
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Eq. (J2Jl is mathematically well-defined on the strip of the 
complex plane determined by the equation < Im(£) < 
(3h, provided that 
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0) < oo, 



(3) 



for all /3i,/?2 > 0. In these conditions, Baym and 
Mermini^ have argued that the function Coit) is ana- 
lytic on the aforementioned domain. In addition, Co(t) 
is uniquely determined on this domain by the values of 
Co{t) on the purely imaginary interval (0, z/3ft), values 



that can be computed efficiently by path-integral Monte 
Carlo techniques (via the Feynman-Kac formula). Fi- 
nally, the correlation function Co(t) is uniquely deter- 
mined at all points of continuity on the frontier of the 
strip < Im(<) < (3h, frontier that obviously includes 
the real axis. 

Berne and Harp 46 have pointed out that the computa- 
tion of thermally-symmetrized quantum correlation func- 
tions 
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with (3 C = (3/2 + it/H and (3 C — (3/2 — it/ 'ft, might be an 
easier computational task, yet the correlation functions 
Go{t) and Co{t) carry essentially the same information, 
because their Fourier transforms satisfy the simple rela- 
tion 

GoH = e-^' 2 C {Lo). 

Certain quantities of physical interest may not even re- 
quire the computation of direct and inverse Fourier trans- 
forms. For the flux (F) autocorrelation functions, Miller, 
Schwartz, and TrompSS have shown that the time in- 
tegrals over the interval [0, oo] of the Cp(t) and Gf(i) 
functions are equal and so, for the determination of the 
quantum rate of reaction)S2i2LS it does not matter which 
of the correlation functions is utilized. As Eq. (J3J) im- 
plies, Go(t) is well-defined in the strip of the complex 
plane defined by the equation |Im(t)| < (3h/2. Baym 
and Mermin's argument can by extended to justify that 
Goit) is analytic in this strip and admits a unique ana- 
lytic continuation from the values of Go (t) on the purely 
imaginary interval {—i(3h/2,i(3h/2). 

Before continuing with our exposition, let us remember 
the statement of the Hamburger moment problem. Sup- 
pose a sequence of real positive numbers {/x^ > 0; k > 1} 
is given. The Hamburger moment problem consists in 
answering the following questions: 

1. Is there a probability distribution dP(tu) on the real 
axis (—oo, oo) such that 



. = / LU k dP{uj), Vfc > I? 



2. If the answer is positive, is the solution unique? (In 
this case, the problem is called determinate.) 

3. If the solution is not unique, can one describe all 
possible solutions having moments 

As a historical note, the moment problem on the interval 
[0, oo) is called a Stieltjes problem, whereas the moment 
problem on a compact interval [a, b] is called the Haus- 
dorff moment problem. The problems are called after the 
names of the mathematicians that have successfully and 
completely resolved the respective problems (the condi- 
tions are slightly different for the three cases, with the 
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Hamburger problem being the most restrictive and chal- 
lenging of the three). For our purposes it suffices to notice 
that a determinate Hamburger solution, if it is a Stielt- 
jes or Hausdorff solution, then it is also determinate in 
the sense of Stieltjes or Hausdorff. Necessary and suf- 
ficient conditions for a sequence of positive numbers to 
be a moment sequence have been given by Hamburger in 
a series of papers from 1920 to 192141 He has also pro- 
duced sufficient and necessary conditions for the problem 
to be determinate. Because the inverse Hamburger mo- 
ment problem lacks continuity with the input data, any 
inversion algorithm designed to recover a probability dis- 
tribution from moment data is ill-conditioned. 

In this section, the computation of thermally- 
symmetrized quantum correlation functions is reduced 
to an inverse Hamburger moment problem. In this re- 
spect, it is first shown that the sequence of derivatives at 
origin of the autocorrelation function is a sequence of mo- 
ments {/Xfc; fc > 0}, up to a normalization factor. In fact, 
the quantum autocorrelation function is the characteris- 
tic function of the probability distribution from which the 
moments fik are derived, probability distribution that is 
commonly called the spectral weight function. The en- 
suing Hamburger moment problem is then shown to be 
determinate. Therefore, the autocorrelation function is 
uniquely determined by its sequence of derivatives at ori- 
gin. While this statement also follows from the Baym and 
Mermin's argument, our proof has the advantage of also 
suggesting reconstruction techniques. As such, Th. [3] 
which identifies a large class of candidate algorithms for 
the inverse Hamburger moment problem, does not follow 
from Baym and Mermin's argument. 



A. The input data 

Let us show that Go(t) is well-defined on the strip 
|Im(f)| < 0H/2 in the complex plane, whenever 

Mo {01,02) = tr {e-^ H OU~^ H 0) < oo, (5) 

for all 0i,02 > 0. With the help of the spectral decom- 
position 



for all 0i,02 > 0. From the inequality 
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Eq. becomes 

G Q {t) = ( [ e-^ E ~^ E '\(E\0\E')\ 2 dEdE', (7) 



whereas the condition given by Eqs. J3J or JSJ now reads 
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one concludes that the integral appearing in Eq. JJJ is 
absolutely convergent for all t > 0, because Go{0) = 
M o {0/2,0/2)<oc. 

In these conditions, Baym and Mermin have argued 
that Go{t) is analytic on the strip |Im(t)| < 0H/2. In 
fact, the analyticity of the autocorrelation function fol- 
lows easily from Eq. JSJ and from the absolute conver- 
gence of the integral appearing in Eq. J7J). Nevertheless, 
for our algorithm, we only need analyticity at origin to- 
gether with a stronger statement on the radius of conver- 
gence of the Taylor series about origin. This is ensured 
by the following proposition. 

Proposition 1 Go{t) is differentiable at origin in- 
finitely many times and the radius of convergence of the 
Taylor series about origin is at least 0H/2. 

Proof. Consider the standard inequality 
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which is valid for all \z\ < r < oo. Pick an arbitrary 
positive number r < K0/2. Then, for all t with \t\ < r, 
we have 
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D k -- 

and 
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The finitude of the last term follows from Eq. JSJ because 
< r < (3h/2. 

Since M r does not depend upon n, an easy inductive 
argument over n and Eq. (|10|) show that the terms Dk 
are finite. Moreover, letting n — > oo in Eq. (|f 01) . we learn 
that 



G Q (t) 



LXJ 



for all i with |t| < r. Since r < /3S./2 is arbitrary, the 
proof is concluded. □ 
From Eq. Q, we notice that G (t) = G Q (-t). There- 
fore, £>2fe+i = for all k > 0. In these conditions, a little 
thought shows that Eq. I|12ll can also be written as 



G {it) 
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fe=0 



(2fc) 



t 2k D 2 k, 



(13) 



the right-hand side series being convergent at least on 
the disc of equation \t\ < h/3/2. Thus, the numbers D 2 k 
are positive [by Eq. l|llfl] and are the even derivatives of 
the imaginary-time correlation function Go (it) ■ 

To summarize, the input data for the algorithm con- 
sidered in the present paper is the sequence of even 
derivatives of the imaginary-time autocorrelation func- 
tion Go (it). This sequence, denoted by D 2 k, consists 
of positive numbers computable by path-integral Monte 
Carlo simulations. 



B. The function that is reconstructed 

The function (distribution) that is reconstructed is the 
power spectrum of the auto-correlation function Goit). 
The power spectrum is defined through the identity 
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and is generally defined as a non-negative tempered dis- 
tribution. With the help of Eq. Q , one computes 
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Simple manipulations lead to 

G (w) = &- M/2 / e- pE \(E + uh\0\E)\ 3 dE, (15) 



which shows that the power spectrum is a non-negative 
distribution. 

By means of Eq. (|14|l . one easily proves that the sym- 
metry of Go(t) implies the symmetry of Go{u). In ad- 
dition, with the help of the inverse Fourier transform 



Go(t) = / e iut G (oj)dw, 



(16) 



one also proves that 
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We summarize the findings of the present subsection 
into the following proposition. 



(12) Proposition 2 The prescription 



dP {u) = —G (u)duj 



(17) 



defines a symmetric probability measure on K. Thus, the 
odd moments fJ-2k+i of the measure are zero. The even 
moments of the probability measure dPo(uj) are finite and 
equal to 

fi 2k = [ L0 2k dP o {u) = Vfc > 1. (18) 

Jr u o 



C. The moment problem to be solved 

Surely, the reader has already anticipated that the 
problem we want to solve is the following Hamburger 
moment problem: Determine the symmetric probability 
measure dPo(uj) on K ; the even moments of which are 
given by the sequence {D2k/Do, k > 1}. However, in or- 
der for the problem to be correctly formulated, we must 
show that there exists a unique symmetric probability 
measure of even moments {Eikl 'Dq, k > 1}. 

The existence is automatically guarantied by the pre- 
scription [Go{u)/ Do]duj, the normalized physical spec- 
tral density, which furnishes an example. For uniqueness, 
we cite the following theorem (Th. 3.11 from Section 2.3 
of Ref.©. 

Theorem 1 // limsup^^^ /2k < oo, then there 
is at most one distribution function Po(w) with fik = 
J oj k dPo(oj) for all positive integers k. 

We then have the following theorem. 

Theorem 2 There exists a unique symmetric probabil- 
ity measure dPo(uj) of even moments {D2k/Do, k > 1}, 
which is the one associated with the physical spectral 
weight function. Consequently, the sequence of positive 
numbers {Z?2fc, k > 0} uniquely determines the autocor- 
relation function Go (t) on the whole real axis. 

Proof Let t = h/3/4 and a = G (it). From Eq. JT^J we 
learn that D 2k < a(2k)\jt 2k . With the help of Stirling's 
formula, we compute 
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(19) 
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and the theorem follows from Th.^and the uniqueness of 
the inverse Fourier transforms of probability distributions 
(so-called characteristic functions of the respective prob- 
ability measures, according to Section 2. 3. a of Ref . |4S|) . 
□ 

In particular, Th. shows that the dynamics on the 
whole line is in principle uniquely determined by the se- 
quence of derivatives at origin of the imaginary-time cor- 
relation function. Of course, this also follows from Baym 
and Mermin's analytic continuation result, but the proof 
we have performed is more direct in the sense that it con- 
nects the uniqueness with the numerical technique in a 
straightforward fashion. The reader will appreciate this 
from the following theorem, which gives general crite- 
ria for the pointwise recovery of the correlation function 
Go(t) on the whole real axis. 

Theorem 3 Let dPo, n (w) be a sequence of symmetric 
probability measures such that 

lim / uj 2k dP .n(u) = D 2k /D 

n — too 

for each k > 1. Then 

lim Gonit) = G (t), yt e K. 

n — >oo 

Observation. Of course, by Go,n(t) we understand, 
up to a multiplication factor of Dq, the characteristic 
function of the measure dPo,n(oj)- The characteristic 
function is defined by 

Go,n(t) =D f e iut dP ,n(")- 

Jg. 

Remembering Eqs. and (|T7Jl . we see that Go(t) is 
also a characteristic function, namely that of the measure 
dPo(uj), because 

Go(t) = D [ e iut dPo(u). 

JR 

Characteristic functions of measures are always continu- 
ous, fact that follows easily from the dominated conver- 
gence theorem. 

Proof of 771.0 Th. 3.12 from Section 2.3 of Ref.|H as- 
serts that the sequence of probability measures dPo t n(<-o) 
converges weakly to dPo(u), because Eq. ifHfl) holds 
true. The first part of the continuity theorem (Th. 3.4 
from Section 2.3 of the same reference) states that the 
weak convergence of the probability measures implies 
pointwise convergence of the corresponding characteris- 
tic functions at all times f e I. The last observation 
concludes the proof of the theorem. □ 

In a sense, Th. [21 says that the pointwise values of the 
correlation functions are the easiest to obtain. Basically, 
any procedure that is capable of reproducing the first 
n moments of the true probability distribution leads to 
convergence of the correlation functions, in the limit of 
large n. Other properties, as for instance certain integral 



values involving correlation functions, are more difficult 
to obtain. Given the general approach put forward in the 
present section, we are now ready to discuss the two main 
computational aspects of the technique: the computation 
of the sequence of even derivatives of the imaginary-time 
correlation function and the numerical resolution of the 
associated Hamburger moment problem. 



III. DERIVATIVES OF THE IMAGINARY-TIME 
CORRELATION FUNCTIONS 

According to Proposition 1, the Taylor series about 
origin of the imaginary-time correlation function Go (it) 
is convergent in the disk of equation \t\ < (3H/2 of the 
complex plane. As the well-known example of the free 
particle flux autocorrelation function (see Ea. l59|) demon- 
strates, in general, one cannot expect convergence be- 
yond this radius. Thus, for the purpose of computing 
derivatives in origin of the imaginary-time correlation 
function, we are forced to restrict the range of values 
of t on which Go(it) is "sampled" to the real interval 
(— Ph/2, (3H/2). On this interval, the correlation func- 
tion Go (it) is computable with the help of the Feynman- 
Kac formulaWSiiS and we now turn our attention to the 
problem of constructing path-integral estimators for the 
evaluation of the high-order derivatives of Go (it) ■ 

We shall illustrate the general strategy for the deriva- 
tion of estimators for the particular case of the flux au- 
tocorrelation function. The reader needs notice that, fol- 
lowing the prescription of Predescu and Doll, 42 we strive 
to bury the time dependence into the potential part of 
the various estimators in order for these estimators to 
have finite variance in the limit of an infinite number of 
path variables. This procedure prevents the well-known 
divergence of the variances of the estimators obtained by 
direct differentiation against imaginary time, with the 
increase of the number of path variables. Such a diver- 
gence is characteristic of the Barker estimators£2i£i and 
is caused by an unfortunate attempt to differentiate the 
Brownian paths entering the Feynman-Kac formula (a fa- 
mous 1933 theorem of Paley, Wiener, and Zygmund says 
that Brownian paths are not differentiable, with proba- 
bility one)!^ In addition, at the cost of utilizing a one- 
dimensional finite-difference scheme, the approach avoids 
the computation of the high order derivatives of the po- 
tential that appear in virial estimators^ as well as in es- 
timators for which the imaginary-time differentiation is 
replaced by the direct action of the Hamiltonian. Even 
more, available numerical results (it is true, for low or- 
der derivatives, only) suggest that the variances of ther- 
modynamic estimators we utilize are smaller than the 
variances for the corresponding virial 43 and Hamiltonian 
techniques^ especially at low temperature. 

For a one-dimensional system, the imaginary-time flux 
autocorrelation function reads^S 

G F (it) = tr ( e -(m+t/h)Hp e -(W-t/h)Hp^ t (2Q) 
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where 



and 



F = [S(x - x s )p+pS (x - x s )] 



2mo 



(21) 



. _ hd_ 

^ i dx 

are self-adjoint operators (therefore, F^ — F). The flux 
operator F corresponds to the dividing surface passing 
through x s (actually, a "dividing point" in this onc- 
dimensional case). Setting fit — 0/2+t/h, Eq. H2U|) takes 
the form 



G F (it) 



( h 



\2rnc 



d 2 p 

p(x,x';P-t)-g^j(x,x';(3 t ) 



2 p 



+ o , (x, x'; (3-t) P (x, x'; fit) 
oxox' 

^X ^'^''^-^ ^ { X ,x';fi t ) 

-^(x,x';0- t )^{x,x';P t ) 



(22) 



X —X — X s 



F 



2^ + R + /3 > 



where, of course, p(x, x'; fit) is the density matrix at the 
inverse temperature fit- 

Let us consider the one-dimensional Feynman-Kac 
formulai^a 



p(x, x'; fi t ) = Pf P (x, x'; fi t )Ee~^ ti v[M«)+**B°]to 

(23) 

which expresses the density matrix as the expected value 
of a functional of the standard Brownian bridge B®. In 
Eq. x r (u) = x+ (x' - x)u and cr t = (h 2 fi t /m ) 1/2 , 

whereas pf p (x,x'; fit) stands for the density matrix of a 
similar free particle at the inverse temperature fit- By 
explicit computation, from Eq. (1221) and the Feynman- 
Kac formula, one derives the equation 



G F {U) = E& ' e -fi-t!^V{x s +a^B<i)du~p t ^ j V{x s +a t Bl')du 



( h 



\2toc 



p fp (0;fi„t)p fp (0;fit)F? (B°,B°') , (24) 



where 



r 1 i r fl 

' V (x s +<T t B°) udu J V (x a + a t B°'^ (1 - u)dv 



+0Lt 



-fi-tfit 
-fi-tfit 



V (x s + <T- t B°) udu 



V' (x s + <r- t B°) (1 - u)du 



J V (x s + a_ t B°)udu J V' (x s + a t B°'}(l-u)du 
V* (x s + <j t B°^ udu J V'(x s +(T-tBl) (l-«)du 



fi-t f V" {x s + a-tB° u )u{l-u)du- fit j V" (x s + a t Bl') u(l - u)du. 



(25) 



In Eq. J23J, the symbols E and E' denote the ex- 
pected values against the independent standard Brown- 
ian bridges B® and B° , respectively. In Eq. l|2*5|l . V'(x) 
and V"(x) denote the first and the second derivatives of 
the potential V(x), respectively. 
Now, Eq. (|24() can be rearranged as 

G F (it) = EE'e- (/3/2) [/o V(x s +a Bl)du+Q V(x s +a Q B°J)du] 



and 



Jo 



V (x s + <7 5„) du 



2 AJ\ {xs + atB o )du _ 



(28) 



Anticipating the use of Monte Carlo techniques for the 
I evaluation of imaginary-time correlation functions and 
T[ (B°,B° ) , (26) related properties, we introduce the normalization factor 

Tfln \ / 



87TTO0 



where 



xe -(0/2)[A_ t (S»)+A t (B«')] 



Mb 



87rmn 



-EE e 



' e -(P/2)[la v { x '+ a oBl)du+^ t V(x s +<j B°J)du] 



(27) 



(29) 

In principle, the factor Mp can be evaluated in a 
separate Monte Carlo simulation, although for the 



8 



one-dimensional example presented later in the paper, 
we shall employ the numerical matrix multiplication 
technique If rate constants rather than absolute 
rates of reaction are desired, one seeks to evaluate the ra- 
tio between J\fp and the partition function of the reactant 
side, Q r . A Monte Carlo approach to the computation 
of such ratios has been recently presented in Ref. 0. 



In any case, the main difficulty in the computation 
of quantum correlation functions does not reside in the 
evaluation of the normalization coefficient Mp- There- 
fore, for the remainder of the present paper, we shall 
focus our attention on the Monte Carlo evaluation of the 
ratios 



G F {it) 



F t (BtB°')) = 



EE e 



> e -(P/2)[Io V(x s +a B°)du+fi V(x s +a Bl')du] ^0 ^ 



/g-C/ 3 / 2 )^ 1 V(x s +a B°)du+^ V(x s +<7„B°')du] 



EE'e 



(30) 



or related quantities. For the purpose of computing av- 
erages of the type given by Eq. I|30|l . it turns out that it 

is useful to replace the estimating function T[ (^B®, B®' 

with the symmetric form 



r t (bi,bi') = \ \tu (j&itf) +n (bib°: 

"(31) 

As follows from the equation Gf(—H) — Gi?(ii), this re- 
placement does not change the value of Gf{H)- However, 
in the next paragraph, we shall prove that the resulting 
estimator has a smaller variance. 

It follows from Eqs. (gSJ and O that 



F'JBl,B°:)=T' t (Bl',Bl 



(32) 



and therefore, 



F t (BlX) = \[H(BtB° 
+Fl(B°,B°')] =f t (B°J,B° 



(33) 



Consequently, the function Tt [B® , B® J is not only sym- 
metric with respect to time inversion, as follows directly 
from Eq. I|3U|) . but also with respect to the exchange of 

variables B® and B®' . Let us write T[ (b®,B®J as the 



sum between its symmetric and its antisymmetric parts 

F t (BlBl')=F t (BlB^) 

'F t (BlBl')-F' t (B^Bi) 



1 

+ 2 



Since antisymmetric functions integrate to zero against 
a symmetric probability measure, and since the products 
of symmetric and antisymmetric functions are antisym- 
metric, we have 



H [BlB* 



= {F t [BlB* 



t'AbIb^-t'Jb^bI 



The last equation and the equality 



T' t B°,B°' = (T t B°,B? = 



G F (U) 
G d (0) : 



which was discussed in the previous paragraph, clearly 
demonstrate that the estimator given by Eq. I|31|) has 
a variance smaller than that of the estimator given by 
Eq. 123). 

To summarize, by Monte Carlo simulations, one may 
compute averages of the type 



G F {it) 
N F 



(b°,b°:)) 



EE / e -(,3/2)[/ 1 V(x 3 +a B°)du+f 1 V(x,+a B°')du] ^ _g0'~j 
EE / e -(/3/2)[/ 1 V(x 3 +<7 B°)du+Q V(x s +a B0')du] 



(34) 



where The estimating function Tt [B®,B®) is symmetric under 



T t (b° 5°') = 1 |jrO f ^ B a B o^ time inversion — that is, T t (#°, = T- t 

"^y/P-tPt — as we u as under the exchange of the variables B® and 



Xe - W2) [A t K)+A_ t ( B r)] + f0 (b°,B°') (35) 

xe -(f)/2)[A- t (B°)+A t (B°')y 
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B®' . Carlo simulations, one may compute the following aver- 

The construction of estimators for derivatives in origin ages 
is straightforward and follows from Eq. (|34|1 ■ By Monte 



1 d k 

7T F 



EE' 



' e -(/3/2)[J^V(x a +aB°)du+S^V(x s +aBZ')dv]^_ :Ft /gO^O'j 



4=0 



t=0 EE'e- {f3/2) y« V(x,+aB°)du+f^ V{x s +aB0J)du] 

I 



(36) 



In this respect, the reader should notice that the function 

T t (B°, B° ,S \ is wcll-dchncd for all t £ (-(3h/2, (3h/2) and 

is infinitely differentiable on this interval provided that 
the potential V(x) is also differentiable infinitely many 
times. In practical applications, the time derivatives ap- 
pearing in Eq. I|36|l are to be computed by finite differ- 
ence. We shall further discuss this matter in Section IV. 

We now describe the construction of estimators for the 
case of a d-dimensional system. For definiteness, we shall 
assume that the physical coordinates have been rescaled 
such that all masses are equal to the common value mo- 
Perhaps after a reorientation of the system of axes so that 
the first coordinate xi is along the reaction coordinate, 
the reactants and products are assumed to be separated 
in the configuration space M. d by a hyperplane of equation 
xi = x s . For the remainder of this section, when dealing 
with expressions involving the density matrix, it turns 
out that it is more convenient to work with the pair of 
position coordinates (x, z), with z = x' — x, rather than 
with the standard (x, x') pair. This is so because identi- 
ties of the type 



dx'pf p (x, x'\ Pt)pf p {x, x'; f3- t )f(x' - x) 
1 



2tt(Tq 



dze- z f(za ±t ), (37) 



^(x,z,5?X' 



~P\ 



where a± t = <7t&-t/coi are clearly simpler to express in 
the new coordinate system. Moreover, transformations 
of the type shown by Eq. Ij37[) are consistent with the 
aforementioned advice of Predescu and Doll that the time 
dependence of paths should be buried into the potential 
part of the Feynman-Kac formula whenever possible. 



With these clarifications, we leave it for the reader 
to demonstrate that the multidimensional analogues of 
the various quantities necessary for the construction of 
derivative estimators are as follows. With the under- 
standing that the quantities V'(x) and V"(x) now denote 
the first order and the second order partial derivatives 
against the reaction coordinate Xi , the multi-dimensional 
analogue of Eq. is 



-p-tpt 
-p-tpt 



J V' (x + <r± t zu + cr t B°'^ udu J V' (x + (T ±t zu + a- t B°'^ (1 - u)du 

[ V (x + a± t z,u + a- t Bl)udu [ V (x + a ±t zu + <T- t B°) (1 - u)d 
Jo J Uo 

J V (x + cr ±t zw + cr_ t B") udu J V' (x + o- ±t zw + a t B^ (1 - u)du 
J V (x + a± t zu + ot-B°') udu J V (x + a± t zu + (T- t B°) (1 - u)du 
P-t J V" (x + a ±t zu + <r- t B°) u(l - u)du ~ Pt J V" (x + a ±t zu + a t B°'^ u(l - u)d 



(38) 



The quantities B® and B®' are independent d- dimensional standard Brownian bridges (ci-dimensional 
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vector valued stochastic processes, the components of The normalization coefficient Afp now reads 
which arc independent one-dimensional standard Brow- 
nian bridges). We also define 



A t (x,z,B°) = f V (x + crozu + cr Q B u )du 
Jo 

[ V (x + o± t zu + a t Bl) du (39) 
Jo 



Wt f 1 



as well as 



, e -(0/2)[A t (x, Z) B°) + A_ t (x, Z ,Br)] + U ZjB l B 0>\ (40) 
Xe -(/3/2)[A_ t (x, Z ,B2)+A t (x, Z ,Bf)] | _ 



J 



Af F = 



1 



1 



87rmo V 27Tf7o 



d-xdzEE' e~ 1 ^ 2 e~ {J3/2) y° V( X +a ^u+a B°)du+J^ V(x+aozu+a B°')du\ 



(41) 



where the integration against the variables x and z is 
done on the (d — 2)-dimensional hyperplane <S, which is 
the subset of the space R d x M. d defined by the equations 
Xt = x s and z\ = 0. Therefore, the symbol dx stands for 
the Lebesgue measure dx2 ■ ■ • dxd, whereas dz stands for 

dz 2 ■ ■ ■ dz d . The Euclidian norm ||z|| = (z 2 H + z^) 1 / 2 

can be replaced by ||z|| = (z| + • • • + z 2 -) 1 / 2 , since the 



coordinate Z\ is kept constant and equal to zero during 
integration. 

In these conditions, up to the value of the normaliza- 
tion coefficient Mf, the derivatives in origin of the flux 
autocorrelation functions can be determined by Monte 
Carlo integration, as implied by the equation 



t=0 



J s dxdzEE'e-INI V (/3/2) [/o V{*+* 0Z u+a a Bl)du+Q V(x+* zu+* B<>')du} |^ / ^ gO^O' 



t=0 



(42) 



IV. SOLVING THE INVERSE MOMENT 
PROBLEM: A NUMERICAL EXAMPLE 

Until now, we have demonstrated that the sequence of 
derivatives at origin completely and uniquely character- 
izes the correlation function. Moreover, the sequence of 
derivatives can be computed by Monte Carlo simulation 
via estimators that have finite variance in the limit of an 
infinite number of path variables (of course, for analytic 
potentials). At this point, it is natural to address the 
problem of recovering the correlation functions from the 
sequence of computed moments. 

More precisely, let us assume that we have com- 



puted the set of even and nonnegative derivatives 
Do, D2, ■ ■ ■ , Dm and that we have calculated the mo- 
ments fX2k — Dzk/Do, for 1 < /c < n. At the very least, 
we would like to construct a sequence of symmetric prob- 
ability distributions dPo,n(u) such that 

M2fe - / 0j 2k dP Otn {uj), (43) 
Jwt 

for all 1 < k < n and n > 1. Indeed, if Eq. (|43|) is sat- 
isfied, then so is the hypothesis of Th. [21 theorem that 
further guaranties that the correlation functions are fully 
recovered (pointwise) in the limit n — > 00. However, 
many times, the pointwise reconstruction of the correla- 
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tion functions does not suffice. For example, in the case 
of the flux autocorrelation function, the chemical physi- 
cists are usually interested in computing the absolute rate 
of reaction, which is the time integral of the correlation 
function 



given by 



k(T)Q r (T) 



G F (t)dt. 



(44) 



Because the first n even moments do not uniquely de- 
termine a symmetric probability distribution, we have 
freedom in choosing the reconstruction algorithm in such 
a way that not only the pointwise values of the correla- 
tion functions, but also various integral expressions are 
recovered in the limit n — > oo. 

Although the optimal reconstruction algorithm de- 
pends upon the nature of the correlation functions and of 
the quantum information being sought, we shall discuss 
and utilize in the present paper a choice that is based on 
the maximum entropy approach. The maximum entropy 
methodi£i22*££i£2*Si suggests that a useful criterion is to 
chose the probability distribution G(u>) that maximizes 
the Shanon entropy 



5(G) 



G(u) In [G(u) /m(u)] du, (45) 



relative to the default model m{u>) and subject to the 
constraints 



G{uj)uj 2k duj = D 2k , 0<k<n. 



(46) 



In information theory, such a probability distribution is 
the least biased one that is compatible with the partial in- 
formation represented by the known first moments. The 
default model m(uj) is a strictly positive distribution. Al- 
though it has a definite probabilistic meaning only if it 
is integrable, non-integrable default models can also be 
used. The choice m(u>) = 1 is called the flat default 
model. 

Simple variational arguments and use of Lagrange mul- 
tipliers show that the unique maximum of the above 
problem is realized for 

Go, n (u) = m(oj) exp ^ XjOJ 2j j . (47) 

The coefficients Ao, . . . , A n are the Lagrange multipliers 
and can be determined from the equations 

D 2k = / m(u))uj 2k exp - V] X 3 uj 2j ] du, < k < n. 

(48) 

Notice that the form of the approximant given by 
Eq. (|47|) ensures both the positivity and the symmetry of 
the power spectrum, properties that have been demon- 
strated in Section II. Then, the entropy of Go,n(w) is 



S [G 



O.n 



Go,n(w) 



x In [Go, n (w)/m(w)] du = ^ X 3 D 2j . (49) 

3=0 

One of the advantages of the maximum entropy algo- 
rithm is that, by use of default models, it may incorpo- 
rate additional physical information that depends upon 
the nature of the quantum results being sought. How- 
ever, for the present example, a flat default model has 
been utilized. Also, for the present application, the data 
have been assumed noiseless. The stability of the final 
results with respect to the errors in the input data has 
been found to be excellent, in part because the number 
of matched moments is small, but also because the dif- 
ferent data are perfectly correlated (they are obtained 
in the same Monte Carlo run). Thus, the assumption 
of noiseless data is good. For larger numbers of included 
moments, more general approaches of Bayesian statistical 
inference with entropic priors also allow for the treatment 
of noise in the data, via likelihood functions^ 

The system of equations (|48f) can be replaced by 



A = In 



— — / m(uj)e 2j=i 
Do 



du 



(50) 



and 



D 2k = A) 



J R m(Lo)cu 2k e-^"^ x ^ 2j du 



1 < k < n. 



(51) 

It is then a simple exercise to verify that Eqs. (|51|) are 
satisfied for all 1 < k < n provided that the Xj 's represent 
the coordinates of the minimum of the entropy functional 



S[G 



O.n 



D In 



Do 



n 
3=1 



(52) 



which is a convex function of Ai, . . . , A„. Due to the con- 
vexity of the function that is minimized, the minimum 
of Eq. (|52[l . if it exists, is unique. The necessary and 
sufficient conditions for the existence of the minimum 
are known in literatureSiSS In the present article, the 
minimization of Eq. Q52JI has been carried out with the 
help of Newton's steepest descent technique. The Hes- 
sian matrix is evaluated explicitly and utilized to predict 
the direction along which to line-minimize. The Golden 
Section search is utilized to optimize along the computed 
direction. As discussed in Ref.|6(J, the computation of the 
coefficients Aj becomes less and less stable as the num- 
ber of matched moments increases and, depending upon 
the number of even derivatives considered, may require 
extended-precision arithmetics. 
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In order to demonstrate its usefulness, we apply the 
moment technique to the problem of computing the quan- 
tum rate of reaction for a symmetric Eckart barrier at 
various temperatures. The parameters for the Eckart 
barrier are chosen to correspond approximately to the 
H + H2 reaction*^ The potential is 



V(x) = Vq sech(ax) 5 



(53) 



with the parameters Vb = 0.425 eV, a — 1.36 a.u., and 
rno = 1060 a.u. 

We evaluate the flux autocorrelation function and its 
first five even derivatives at origin by Monte Carlo simu- 
lations, as described in Section III. The derivatives of the 
estimator T t (B®, B®) appearing in Eq. I|36|l are replaced 
by numerical approximations computed via central differ- 
ence. Remembering that Tt{B®, B®') is symmetric under 
the transformation 1 1— ► —t, the finite-difference formulas 
take on the general form 



dt 2k 



3=0 



(54) 

where the coefficients Cj t % are given in TableQ] Numerical 
experiments demonstrate that a time step of 



64~T 



is sufficient for a determination of the derivatives to an 
accuracy of less than 2%. 

Regarding the computation of derivatives by finite dif- 
ference, the range of values of t that can be utilized 
depends on the order of the derivatives as well as on 
the numerical precision with which the computations are 
conducted. For the present paper, we employ the IEEE 
floating-point data type double (64 bit) for the repre- 
sentation of real numbers. Increasing the order of the 
derivatives beyond 10 requires use of extended-precision 
data types' 



2Jb 


Cfe,0 


Ck.l 


Ck,2 


Cfc,3 


Cfc,4 


Cfc.5 





1 

















2 


-5269/1800 


10/3 


-10/21 


5/63 


-5/504 


1/1575 


4 


1529/120 


-1669/90 


4369/630 


-541/420 


1261/7560 


-41/3780 


6 


-1023/20 


323/4 


-39 


87/8 


-19/12 


13/120 


8 


154 


-252 


136 


-46 


26/3 


-2/3 


10 


-252 


420 


-240 


90 


-20 


2 



TABLE I: Numerical values for the coefficients Chj appearing 
in the finite-difference approximations of the derivatives of 
order 2k. 

Many times, the chemical physicist takes the differ- 
ent approach of constructing models (and, therefore, em- 
pirical inversion techniques) that have already incorpo- 
rated additional physical inputS In such cases, the fi- 
nite number of derivatives that can be computed using 
the data type double may suffice for many practical pur- 
poses. This is why it is appropriate to table the coeffi- 
cients Cj t k, for the reader's convenience. General rules for 



computing derivatives of arbitrary orders and accuracy 
have been discussed elsewhere*^ According to Eq. 
the accuracy of the finite-difference scheme is largest for 
the small-order derivatives and decreases for the larger- 
order derivatives, if all the information contained in the 
6 points at which Tt{B®, B® ) is evaluated is to be taken 
into consideration. This is to our advantage, because the 
low-order derivatives are computed with increased preci- 
sion despite the relatively large value of the discretization 
step r demanded by the higher-order derivatives. 

For the sake of an example, in Table II, we present the 
Monte Carlo estimates of the first five even derivatives at 
origin for the Eckart barrier at the temperature of 100 K. 
The derivatives have been evaluated in 10 million Monte 
Carlo points with the help of the estimators introduced 
in Section III. For the discretization of the Feynman-Kac 
formula, we employ Predescu's fourth-order path-integral 
technique^ 5 , with a number of 64 path variables. This 
technique is basically a Trotter product 



p n (x,x';/3) = [ dxi . . . [ dx n po (x, xi; — 

Jr Jr V n + 1 



' f 3 

. po I x n ,x ; 



71 + 1 



(56) 



(55) of a short-time approximation of the type 



p (x,x';(3) = pf p (x,x';0) / dp,(ai) ■ ■ ■ / dp(a q ) 



exp < -0^2 w i V 
{ »=i 



k=l 



(57) 



The quadrature points Ui and weights Wi as well as the 
functions Afc(it) are designed such that the convergence 

p n (x,x';f3) -> p{x,x';/3) 

is as fast as 0(l/n 4 ). These parameters are universal, 
in the sense that they are independent of the choice of 
potential V(x), and are given in Ref. l65l reference that 
should be consulted for further information. 

At this low temperature of 100 K, the Monte Carlo 
sampling requires the use of parallel temperingj&S 7 . 
which, however, successfully copes with the sparse sam- 
pling problem caused by the crossing and recrossing of 
the barrier by the Brownian paths. As a matter of 
fact, by Monte Carlo integration, we compute the ra- 
tios D2k/J^F and the associated statistical errors (two 
standard deviations). The quantity j\fp is evaluated 
with the help of the numerical matrix multiplication 
technique which provides essentially exact results. 
Thus, the relative errors reported in Table II are equal 
to the relative errors of the ratios D^k/J^F and are, there- 
fore, representative of the variances of the estimating 
functions utilized in the Monte Carlo simulation. 

Once the power spectrum Gf,u{u) is determined, the 
absolute rate of reaction can be computed from Eq. 144fl , 
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Order 





2 


4 


6 


8 


10 


Value 


5.787E-17 


2.389E-22 


4.010E-27 


1.395E-31 


7.985E-36 


6.781E-40 


Error 


2.5% 


2.4% 


2.4% 


2.7% 


3.9% 


6.1% 



TABLE II: Derivatives (second row) and relative errors (third 
row) for the symmetric Eckart barrier at 100 K. The errors 
are twice the percentile relative value of the standard devia- 
tion. The errors do not include the systematic errors due to 
the utilization of finite-difference approximations, which have 
been estimated to increase the final errors with less than 2%. 



as the quantity 

k(T)Q r (T) = 

_ 1 
~ 2 

Let us remember that 
G Fn {t) - 



G Fn (t)dt 



G F , n (t)dt = wG F , n (p). (58) 



G F (t), vt e 



for all reconstruction algorithms that satisfy the hypothe- 
sis of Th.|31 However, as already mentioned several times, 
this does not automatically imply pointwise convergence 
in the frequency domain. Sure enough, convergence in 
the frequency domain is necessary only for the purpose 
of computing the absolute rate of reaction as the time 
integral of the flux autocorrelation function, the power 
spectrum of which is continuous at origin. It is not re- 
quired for other autocorrelation functions. Because it 
depends on the physical significance of the correspond- 
ing autocorrelation functions and on the nature of the 
quantum information that is sought, the development of 
optimal reconstruction algorithms is a case by case prob- 
lem. 

It is beyond the purpose of this paper to conduct any 
mathematical proofs related to the pointwise convergence 
of the power spectrum of the flux autocorrelation func- 
tions. However, the percentile relative errors for the ab- 
solute rates of reaction presented in Table III strongly 
suggest that the maximum entropy algorithm discussed 
in previous paragraphs is viable for the purpose of com- 
puting rates of reaction. The errors eventually increase 
as the temperature is lowered, but the reader may notice 
that the relative errors are sufficiently small to make the 
algorithm useful even in the tunneling regime of temper- 
atures (T < 300 K). 

At large temperatures, the relative errors converge to 
the relative errors for a free particle. The thermally- 
symmetrized flux autocorrelation function for the free 
particle io 22 i 38 



G F {t) 
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Its power spectrum reads 
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(3h 2ir 



Luhf3 



(59) 



(60) 



where K\ (x) denotes the respective modified Bessel func- 
tion of the second kind. The function xKi(x) is contin- 
uous at origin, indeed, but its even derivatives in origin 
are not defined. Therefore, the function xK\(x) is not 
readily approximated around origin by smooth functions 
of the type given by Eq. (|47p. Thus, for example, a use- 
ful direction for future research is to modify the default 
model in the the maximum entropy algorithm so that to 
properly account for the known high-temperature limit. 



Order of 
derivatives 


Temperature 


100 K 


200 K 


300 K 


500 K 


1000 K 


2000 K 


oo 


2 


-13.8 


-2.3 


8.4 


-2.1 


-18.3 


-25.7 


-27.6 


6 


-4.9 


-0.8 


2.5 


1.8 


-7.7 


-15.0 


-17.1 


10 


-2.9 


0.3 


0.0 


1.3 


-5.4 


-11.9 


-13.4 



TABLE III: Percentile relative errors for the absolute rates 
of reaction computed using all derivatives up to the maxi- 
mum orders of 2, 6, and 10, respectively. The errors are given 
as functions of temperature. Whenever the minimization al- 
gorithm did not converge properly while using the maximal 
number of derivatives, a smaller number of derivatives has 
been utilized. The relative errors for the high-temperature 
limit are those for the free particle case (which are indepen- 
dent of temperature). 



V. SUMMARY AND DISCUSSION 

A new technique for extracting quantum dynamical in- 
formation from imaginary-time data has been proposed. 
The technique consists in solving a symmetric Hamburger 
moment problem with even-order moments related to the 
even-order derivatives at origin of the quantum autocor- 
relation function. It has been demonstrated that the 
derivatives at origin uniquely determine the autocorre- 
lation function. The derivatives can be computed by 
Monte Carlo simulations with the help of estimators of 
finite variance. The pointwise reconstruction of the auto- 
correlation functions can be performed by those inversion 
algorithms that satisfy the hypothesis of Th.|3 although 
additional care is needed if other quantities, as for in- 
stance certain integral values, are also sought. A moment 
based maximum entropy inversion algorithm has been 
numerically shown to cope successfully with the problem 
of computing absolute rates of reaction for a symmetric 
Eckart barrier. 

Perhaps, the most important step in the present de- 
velopment is the realization that the derivatives at origin 
of the imaginary-time autocorrelation functions are com- 
putable solely by Monte Carlo simulations. As argued 
in the introduction, the sequence of derivatives at origin 
represents a set of data that is more suitable for the prob- 
lem of extracting quantum dynamical information than 
the mere Monte Carlo evaluation of the imaginary-time 
autocorrelation function on a grid. However, future re- 
search is necessary in order to quantify in precise manner 
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the efficiency of the new algorithm. In particular, the 
scaling of the variances of the Monte Carlo estimators 
with the degree of the derivatives, the dimensionality of 
the physical system, and the temperature must be deter- 
mined. 

The numerical results presented in Section IV demon- 
strate that the derivatives at origin of autocorrelation 
functions contain useful information that can be utilized 
in at least two ways. First, one may employ this in- 
formation together with various inversion algorithms for 
the Hamburger moment problem. In this respect, I be- 
lieve that methods of Bayesian statistical inference and 
maximum entropy will be most useful, especially because 
such techniques can incorporate additional physical infor- 
mation (as, for instance, a certain limiting behavior) by 
appropriate choices of default models. Second, if only a 
small number of derivatives are computed, the chemical 



physicist has also the option of developing certain phys- 
ical models depending on parameters that can be deter- 
mined from matching the known derivatives. Which of 
these two ways will be the most successful for practical 
applications remains to be seen. 



Acknowledgments 

The author acknowledges support from National Sci- 
ence Foundation through Grant No. CHE-0096576. He 
wishes to express a special thanks to Professor William 
H. Miller for suggestions and stimulating discussions 
concerning the present development. He also acknowl- 
edges an anonymous referee, whose useful comments have 
helped improve the presentation of the paper. 



* Electronic address: cpredescu@comcast.net 

1 R. P. Feynman, Rev. Mod. Phys. 20, 367 (1948). 

2 A. M. Amini and M. F. Herman, J. Chem. Phys. 99, 5087 

(1993) . 

3 B. J. Berne and D. Thirumalai, Annu. Rev. Phys. Chem. 
47, 401 (1986). 

4 R. H. Cameron, J. Math. Phys. 39, 126 (1960). 

5 J. D. Doll, J. Chem. Phys. 81, 3536 (1984). 

6 V. S. Filinov, Nucl. Phys. B 271, 717 (1986). 

7 N. Makri and W. H. Miller, Chem. Phys. Lett. 139, 10 
(1987). 

8 J. D. Doll, T. L. Beck, and D. L. Freeman, J. Chem. Phys. 
89, 5753 (1988). 

9 C. H. Mak, Phys. Rev. Lett. 68, 899 (1992). 

10 G. Baym and D. Mermin, J. Math. Phys. 2, 232 (1961). 

11 E. Nelson, J. Math. Phys. 5, 332 (1964). 

12 H.-B. Schiittler and D. J. Scalapino, Phys. Rev. Lett. 55, 
1204 (1985). 

13 S. R. White, D. J. Scalapino, R. L. Sugar, and N. E. Bick- 
ers, Phys. Rev. Lett. 63, 1523 (1989). 

14 M. Jarrell and O. Biham, Phys. Rev. Lett. 63, 2504 (1989). 

15 M. Jarrell and J. E. Gubernatis, Phys. Rep. 269, 133 
(1996). 

16 J. E. Gubernatis, M. Jarrell, R. N. Silver, and D. S. Sivia, 
Phys. Rev. B 44, 6011 (1991). 

17 D. Thirumalai and B. J. Berne, J. Chem. Phys. 79, 5029 
(1983). 

18 E. Gallicchio and B. J. Berne, J. Chem. Phys. 101, 9909 

(1994) . 

19 D. Kim, J. D. Doll, and J. E. Gubernatis, J. Chem. Phys. 
106, 1641 (1997). 

20 D. Kim, J. D. Doll, and D. L. Freeman, J. Chem. Phys. 
108, 3871 (1998). 

21 G. Krilov, E. Sim, and B. J. Berne, J. Chem. Phys. 114, 
1075 (2001). 

22 W. H. Miller, S. D. Schwartz, and J. W. Tromp, J. Chem. 
Phys. 79, 4889 (1983). 

23 E. Rabani, G. Krilov, and B. J. Berne, J. Chem. Phys. 
112, 2605 (2000). 

24 E. Sim, G. Krilov, B. J. Berne, J. Phys. Chem. A 105, 
2824 (2001). 



25 E. Rabani, D. R. Reichman, G. Krilov, B. J. Berne, P. 
Natl. Acad. Sci. USA 99, 1129 (2002). 

26 S. R. White, Phys. Rev. B 44, 4670 (1991). 

27 S. R. White, Phys. Rev. B 46, 5678 (1992). 

28 M. Caffarel and D. M. Ceperley, J. Chem. Phys. 97, 8415 

(1992) . 

29 J. Deisz, K.-H. Luk, M. Jarrell, and D. L. Cox, Phys. Rev. 
B 46, 3410 (1992). 

30 J. Deisz, M. Jarrell, and D. L. Cox, Phys. Rev. B 48, 10227 

(1993) . 

31 J. Skilling, in Maximum Entropy and Bayesian methods, 
edited by J. Skilling (Kluwer, Dordrecht, 1989), p. 455. 

32 R. N. Silver and H. Roder, Int. J. of Mod. Phys. C 5, 735 

(1994) . 

33 D. A. Drabold and O. F. Sankey, Phys. Rev. Lett 70, 3631 
(1993). 

34 L. W. Wang, Phys. Rev. B 49, 10154 (1994). 

35 R. N. Silver, H. Roder, A. F. Voter, and J. D. Kress, J. of 
Comput. Phys. 124, 115 (1996). 

36 R. N. Silver and H. Roder, Phys. Rev. E 56 4822 (1997). 

37 W. H. Miller, J. Chem. Phys. 61, 1823 (1974). 

38 W. H. Miller, J. Phys. Chem. A 102, 793 (1998). 

39 G. A. Athanassoulis and P. N. Gavriliadis, Prob. Engng. 
Mech. 17, 273 (2002). 

40 E. T. Jaynes, Phys. Rev. 106, 620 (1957). 

41 A. Tagliani, J. Math. Phys. 35, 5087 (1994). 

42 C. Predescu and J. D. Doll, J. Chem. Phys. 117, 7448 
(2002). 

43 C. Predescu, D. Sabo, J. D. Doll, and D. L. Freeman, J. 
Chem. Phys. 119, 12119 (2003). 

44 J. N. Lyness, Math. Comput. 22, 352 (1968). 

45 J. D. Doll, M. Eleftheriou, S. A. Corcelli, and David L. 
Freeman, Quantum Monte Carlo Methods in Physics and 
Chemistry, edited by MP. Nightingale and C.J. Umrigar, 
NATO ASI Series, Series C Mathematical and Physical 
Sciences, Vol. X, (Kluwer, Dordrecht, 1999). 

46 B. J. Berne and G. D. Harp, Adv. in Chem. Phys. 17, 63 
(1970). 

47 H. Hamburger, Math. Ann. 81, 235 (1920); 82, 120 (1921); 
82, 168 (1921). 

48 R. Durrett, Probability: Theory and Examples, 2nd ed. 



15 



(Duxbury, New York, 1996). 

49 B. Simon, Functional Integration and Quantum Physics 
(Academic, London, 1979). 

50 J. Barker, J. Chem. Phys. 70, 2914 (1979). 

51 M. F. Herman, E. J. Bruskin, and B. J. Berne, J. Chem. 
Phys. 76, 5150 (1982). 

52 R. Paley, N. Wiener, and A. Zygmund, Math. Z. 37, 647 
(1933). 

53 C. Predescu, D. Sabo, J. D. Doll, and D. L. Freeman, J. 
Chem. Phys. 119, 10475 (2003). 

54 A. D. Klemm and R. G. Storer, Aust. J. Phys. 26, 43 
(1973). 

55 D. Thirumalai, E. J. Bruskin, and B. J. Berne, J. Chem. 
Phys. 79, 5063 (1983). 

56 T. Yamamoto and W. H. Miller, J. Chem. Phys. 120, 3086 
(2004). 

57 E. T. Jaynes, in The Maximum Entropy formalism, edited 
by R. D. Levine and M. Tribus (MIT Press, Cambridge, 
1978), pp. 15-118. 

58 J. Skilling, in Maximum Entropy and Bayesian methods, 
edited by J. Skilling (Kluwer, Dordrecht, 1989) p. 45. 



59 S. F. Gull, in Maximum Entropy and Bayesian methods, 
edited by J. Skilling (Kluwer, Dordrecht, 1989) p. 53. 

60 A. Tagliani, J. Comput. Appl. Math. 90, 157 (1998). 

61 W. H. Miller, Y. Zhao, M. Ceotto, and S. Yang, J. Chem. 
Phys. 119, 1329 (2003). 

62 Y. Hida, X. S. Li, and D. H. Bailey, in Proceedings of the 
15th IEEE Symposium on Computer Arithmetic, IEEE 
Computer Society, 2001, pp. 155-162. 

63 N. F. Hansen and H. C. Andersen, J. Chem. Phys. 101, 
6032 (1994). 

64 I. R. Khan and R. Ohba, J. Comput. Appl. Math. 107, 
179 (1999). 

65 C. Predescu, Phys. Rev. E 69, 056701 (2004). 

66 C. J. Geyer, in Computing Science and Statistics: Proceed- 
ings of the 23rd Symposium on the Interface, edited by E. 
M. Keramigas, (Interface Foundation: Fairfax, 1991), pp. 
156 - 163. 

67 K. Hukushima and K. Nemoto, J. Phys. Soc. Jpn. 65, 1604 
(1996). 



